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Abstract 



The Collective Beam-Beam interaction is studied in the framework of maps with a 
"kick-lattice" model in the 4-D phase space of the transverse motion. A novel approach to 
the classical method of averaging is used to derive an approximate map which is equivalent 
to a flow within the averaging approximation. The flow equation is a continuous-time 
Vlasov equation which we call the averaged Vlasov equation, the new model of this paper. 
The power of this approach is evidenced by the fact that the averaged Vlasov equation has 
exact equilibria and the associated linearized equations have uncoupled azimuthal Fourier 
modes. The equation for the Fourier modes leads to a Fredholm integral equation of the 
third kind and the setting is ready-made for the development of a weakly nonlinear theory 
to study the coupling of the tt and a modes. The it and a eigenmodes are calculated from 
the third kind integral equation. These results are compared with the kick-lattice model 
using our weighted macroparticle tracking code and a newly developed, density tracking, 
parallel, Perron-Frobenius code. 

PACS: 02.30.Rz, 29.20.Dh, 29.27.Bd, 45.20.Jj, 52.59.-f, 52.65.-y 
Submitted to New Journal of Physics 

1 Introduction 

In this paper we introduce a new model for the collective beam-beam interaction for hadron 
beams in 4D transverse phase space (2 degrees-of-freedom) . This both generalizes and simplifies 
the work of [H EJ [3j H] on the collective beam-beam interaction in high energy colliders. In 
addition, it extends the preliminary 1 degree-of- freedom collective case in [5]. Our model is 
based on the classical method of averaging generalized to maps and collective forces. We do not 
distribute the beam-beam force around the ring as is usually done. The technique we introduce 
should be of general interest for studies of Vlasov systems with a localized perturbative collective 
force. 

In Section EJ we discuss our basic kick-lattice model for the evolution of the 4D phase space 
densities of the two beams. In Section [31 we briefly review the basic averaging theory which 



is generalized in this paper. Previously, this averaging theory was applied to the weak-strong 
beam-beam in one and two degrees-of- freedom [HI E] . The equations of the kick-lattice model 
will be transformed to a standard form for the method of averaging in Section H] and the general 
"averaged Vlasov equation" (AVE) will be derived (See equation fl2~4l) ). We then introduce the 
special case we treat in this paper, namely the case where the tunes of the two beams are 
identical and are non-resonant. In this case the AVE has the property that any function of 
the action only is an equilibrium solution. In Section we linearize about these equilibria and 
discuss the linearized equations and the associated third kind integral equation. In addition, 
we compare our integral equation with the analogous integral equations which arise in the 
standard plasma problem and in the beam dynamics problems concerning the longitudinal 
dynamics with wake fields for a coasting beam and a bunched beam. In Section [6] we present 
numerical results for the it and a mode eigen-problems for an axially symmetric Gaussian 
equilibrium and compare these results with simulations on the exact model of Section [2J In 
Section 7 we give a summary and point to future work. An appendix is included which gives a 
first principles calculation of the beam-beam force. 



2 The Kick-Lattice Model in 4D Phase Space 

To describe the evolution equations, we refer to the bunches as "unstarred" and "starred" , and 
for every quantity X describing the unstarred bunch, the quantity X* describes the starred 
bunch. The evolution equations are symmetric: the equation for the starred bunch is obtained 
from the unstarred bunch by interchanging starred and unstarred quantities, so we mostly state 
only the equation for the unstarred bunch. 

We consider two counter-rotating particle bunches, which collide head on at a single interac- 
tion point (IP). The electromagnetic interaction at the IP is determined up to a proportionality 
factor by the dimensionless "potential" $, which satisfies the Poisson equation — A$ = 2ixp* . 
Here p*(x,y) is the spatial density (normalized to one), and the potential $[p*] : M. 2 — > R is 
given by 

®[p*](x,y) = / G{x-x',y-y')p*(x',y')dx'dy', (1) 



where G(x,y) = — hi(y/x 2 + y 2 /cr) = — ln(>/x 2 + y 2 ) + ln(cr) is the Green's function. In 
the following we will omit the scaling factor o which is in principle needed for dimensional 
correctness but which can be chosen completely arbitrarily since it does not contribute to the 
beam-beam kick. 

Letting n refer to the state of the system just before the IP, particles in the unstarred bunch 
change their phase-space position u = (x,y,p x ,p y ) T according to the map 

U n+ 1=M ^n + C^O,^^]^),^^]^)^ ). (2) 

The associated phase space density ip n evolves via ip n+ i(u n+ i) = ip n (u n ), c 

d „ , d x ' ' 



Mu) = ^n+i | M \ u + C ( 0,0,— $[p* n )(Vu),— $[p* n ](Vu) ) ) ) . (3) 



which is easily inverted to give 



C (o, 0, ^[p^VM^u), ^[^(VM-'u^j j . (4) 
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Here Ai is a stable linear symplectic map representing the linear lattice, ( is the beam-beam 
factor, V = ^ * ° ° ° J projects phase space on configuration space, and the spatial and 

phase-space densities are related by 

P*(x,y)= 4>*(x,y,p x ,p v )dp x dp y . (5) 

The beam-beam factor, which is derived in the appendix, is ( = ^~ r p; where the absolute 

value of r p = 47r ^ me2 is the classical particle radius (as long as only elementary particles or ions 
of the same charge state are involved), N is the number of particles, q is the particle charge, 7 
is the Lorentz factor associated with (3, and m is the particle mass. For all modern colliders, i.e. 
in the limit (3,(3* — > 1, £ can be approximated by £ rs ^—r p . The evolution law for the starred 
beam is obtained by replacing .M by by ^ and £ by C* where starred and unstarred are 

interchanged in (3, 7, N and m. 

Equation d2J) can be written more compactly as 



u n +i 



M(u n + CJ^ u ®[ P * n ]{Vu n )), (6) 



where J^fc = ( _° fe j is the unit symplectic matrix. We note here that a map is said to be 

symplectic if the Jacobian, M, of the map satisfies M T JM = J . We have written the kick in 
a "Hamiltonian form" because eventually a transformed $ will be a Hamiltonian for a flow. 
For simplicity, we take 
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where Cj := cos(27rz/j), S*j := sin(27r^j) and where z/j for i — x,y are the tunes. We have assumed 
that the beta functions, (3 X and (3 y , have minima at the IP. The distinction between the lattice 
(3 and the relativistic (3 should be clear from context. 

To relate ( to the usual beam-beam parameter, we linearize the kick in (j2J) about (x, y) = 
(0, 0) in the case where p* is mirror symmetric and invariant under | rotations, i.e. when 
p*(x, y) = p*(—x, y) = p*(x, —y) = p*(y, x). Note that this is still a weaker constraint than full 
axial symmetry. Because of these symmetries, $ z (0) = & y (Q) = $xy(0) = and $ xa; (0,0) = 
$ OT (0, 0) = — 7rp*(0, 0) where the latter uses Poisson's equation. Thus the kick matrix becomes 
f *j ° T \ where k = 7rCp*(0,0). The tune shift is & = = —^Pik + 0{k 2 ) with i = x,y. 

Thus the beam-beam parameter & = — j(3i(p*(0, 0). For a round Gaussian, this gives the 
standard result. 



3 Map Averaging and Error Bounds 

Here we give an overview of the averaging formalism, which we generalize in this paper, and 
briefly discuss error bounds. We consider the autonomous "kick-rotate" map in M 2 

u n+l = (u n + e(0, -$'Kn)) T ) (8) 
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with the small parameter e. This is a model for the one degree of freedom weak-strong beam- 
beam interaction and was discussed in [5J. The transformation 



u = 



V 



(9) 



leads to the non-autonomous map 



v n +i = v n + eJ 2 V v H(v n , nu), 



(10) 



where H(v,8) = V(v\ cos(2tt$) + t> 2 sin(27r#)). This is in a standard form for the method of 
averaging in which the transformed dependent variable, v, is slowly varying. If v is irrational, 
then from Weyl's equidistribution theorem [Hj the average of H(v, nv) over n exists and is given 
by H(v) = Jg 1 H(v, 6) d6. It is therefore natural to ask, for what values of v are solutions of 
(ITUi) approximated by solutions of the averaged map 



Even though the maps in (jBJ) and (TTUT) are symplectic, the averaged map is not. However the 
averaged flow associated with (fTTl) and defined by 



is Hamiltonian, and it is easy to show that \w n — w(n)\ = 0(e) over 0(1/ e) times. Since ( fl2|) is 
autonomous, (TTTj) can be viewed as the Euler method for numerically integrating ([121) . Approx- 
imating equation (flOj) with (fTTl) is considered in our previous work [6J, where we introduce the 
concept of a far-from-low- order-resonance zone for v. This zone is formed by removing a finite 
number of intervals centered on low-order rationals, therefore v needs to satisfy only finitely 
many Diophantine conditions, and does not need to be irrational, which makes the formalism 
much more useful in the applications. The error bound \v n — w n \ = 0(e) is obtained without the 
usual near-identity-transformation and is unchanged asymptotically if an 0(e 2 ) term is added 
to equations (f8]fT0l) . 

In [7] , we extend the formalism of equations (l8l fT2l) to the weak-strong beam-beam interac- 
tion in 2 degrees-of-freedom. The equation of motion corresponding to (JED is just equation (J2j) 
with p* replaced by the spatial density of the strong beam. We are working out the details of 
the averaging theorem in this more complicated case with two frequencies [9]. Some ingredients 
of our approach can be found in [TDj. We generalize this to the collective beam-beam interaction 
in the next section. 

4 Map- Averaging for Vlasov Systems 

We will begin by transforming (j2J) using a representation of the solution to the unperturbed, 
C = 0, problem. The new coordinates will be slowly varying if ( is small. As in the previous 
section we could proceed by letting u = A4 n q which gives 



w n+1 = w n + eJ 2 V w H(w n ). 



(11) 



w = eJ 2 V w H(w) 



(12) 



q n +i = q n + (M n J 4 



r iV u ^[p* n ](VM n q n ). 



(13) 



The averaged equation then becomes 



w n+1 = w n + QJaV w F[p*\(w n ) 



(14) 
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where F[p*}(w) denotes the n-average of M~ n J^V u ^[p*}(VM n w). 

The transformation to (1131) turned out to be a major advance in Sobol's implementation of 
the Perron- Frobenius (PF) method [TT] (See [12] and [13] for a discussion of the PF method). 
In addition, ffl3l) is well suited for an error analysis which is in progress [Hj . However, an 
action-angle transformation may be better suited to understand approximate equilibria and 
the associated linear analysis that we do here and that is how we will proceed. 

The action-angle transformation from u = {x,y,p x ,p y ) T to slowly varying coordinates v = 
(Q x , Q y , J x , J y ) T is given by 



x 

Px 



\/2J x p x sm(27rnu x + Q x ) (15) 
^2J X /(3 X cos(2imv x + 6^) (16) 



y = J 2 J y j3 y sm(27rnp x + Q y ) (17) 



Py = j 2Jy/ j3yCOS(2ltni>y + Q y) . (18) 

Note that for fixed J and 6 these are solutions of the equations of motion with ( = 0, that is 
without the beam-beam force. 
Equation ([2]) becomes 



where 



Vn+i = v n + CJiVvH[**„\(v n ,n) + 0(C), (19) 

H[V*]{v,n) := / $*(v')dv' 

x G^y/2p x J x sm(2imv x + Q x ) - a/2/5* J' x sin(27mz/* + 6^.), 

y/2/3 y J y sm{2Trnu y + Q y ) - ^2(3* J' y sin(27rn^ + ej,)). (20) 

The integral in ( 1201) is taken over [0, 2ir] in the G's and over [0, oo) in the J's. Since the 
transformation is symplectic it is also volume preserving. Thus the (G, J)-density is given by 
*&n( v ) = 4>n(u), and its evolution law is 

= V n+1 (v + (J 4 V v H[$* n }(v,n) +0(C 2 )), (21) 

or equivalently 

V n+1 (v) = V n (y - (J 4 V v H[V* n ](v,n) +0(C 2 )). (22) 

Clearly v n and v* are slowly varying for ( and small, and it follows that the transformed 
densities \& and are slowly varying. Thus ffl9|) is in a standard form for averaging and 
we now follow the procedure laid out in the previous section. The averaged map problem is 
obtained from (TEfl) by replacing H by the appropriate n-average H and dropping the 0(( 2 ) 
term. The associated averaged flow problem is autonomous and has the Hamiltonian form 

w = (J 4 V w H[**]{w). (23) 

Thus the averaged Vlasov equations for \l/ and \&* become 

$¥ + C{¥, = 
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where {f,g} = ^f^r + £t4t ~ If Ir - Iflr is the Poisson bracket. Note that H*W] is 

U, ' ,J O0X OJx O0y OJy dJx oOx dJy Otiy L J 

obtained from (|20|) by interchanging the starred and unstarred parameters (3 and v. System 
(1240 is the new model referred to in the title. Since £ and are small one immediate advantage 
of (|24p over §Z§ is that the step size in a numerical integration of (1241) can be 0(l/max(£, £*)) 
which is much larger than one turn. 

At this stage the problem is general with parameters (v x , u y , u*, u*, f3 x , /3 y , j3 x , /3* (, (*) and 
the correct averaged Hamiltonian H depends on the relation between the four tunes. Here we 
discuss the case v x = u* and u y = v* because (i) we wish to compare and contrast our results 
with 0, H] and (ii) it simplifies the calculation of the average. In this case G(- ■ •) in ( |20l) can 
be rewritten as G(D X sm(2irnu x + (p x ), D y sin(2imUy + ip y )), where 

D x = D(p x j x ,p:,j' x ,e x -e' x ), 

D y = D( f3 y J y ,(3;j' y ,e y -Q' y ) , (25) 
D(r, s, t) = \j2r + 2s - 4^ cost, (26) 

and the phases <p x and (p y are easily determined from the trigonometry involved. If in addition 
we consider the case where v x and v y are non-resonant (in the sense that k x v x + k y v y = k Q 
k x = k y = ko = 0), then the averaging over nv x and nv y can be done separately and each 
average can be replaced by the associated integral. Thus the averaged Hamiltonian becomes 

H[**](v)= [ G(D x ,D y W(v')dv', (27) 

JT 2 xR 2 + v j 

where G(D x ,D y ) := l/(27r) 2 L 2 G(D x smt x , D y sint y ) dt x dt y and as before v = (Q,J). From 
our experience with the non- collective case [6J EJ, we expect this to be valid for u x , v y far from 
low-order resonances; work on the error estimates is in progress [2]. Note that H = H* if 

(3 X = (3* and f3 y = f3* y . 

Because of the convolution structure of the 9 integral in (J27J) (see (125]) ) any function \1>* = 
\I/*(J) results in H being independent of 0. It follows that any pair of densities \l/ e and 
that are independent of 9 are an equilibrium pair for (1241 . Since ipn(x,y,p x ,p y ) = \l/n(@, J) ~ 

\& e (J) = \& e ^| (j^ + P x p x ^j , \ (j^ + PyPl^, we see that (T5]) and equivalently ( jl~9{) have quasi- 

equilibria given the averaging approximation. In [5], we have verified this in the 2D phase 
space case. Figure [T] shows the evolution of two action densities at a hundred out of more 
than 10 5 turns under the exact map, the 2D analog of (j4]), using the PF method for tracking 
phase space densities. The parameters are v — v* = VE - 2 and f = 3 x 10~ 3 . The red 
crosses represent the evolution of a quasi-equilibrium, namely the centered Gaussian \E'o(@ ) J) — 
^e(J) = l/(2-KE)e~ J / £ (for e = 1), while the green x's stand for a Gaussian that was initially 
shifted by la x , giving a 6-dependent density. Since the red crosses for different discrete J lie on 
top of each other, \l/ e hardly evolves over the 10 5 turns which is consistent with the averaging 
result. The O-dependent density however, fluctuates as is indicated by the band of green x's. 
Thus we have strong evidence for the existence of quasi-equilibria, however we believe that (j3j) 
does not admit exact equilibria. This is in contrast to the lepton case, see [T5] . 
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p n (J): CR / ^=.003 / AQ=0 / PF 
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Figure 1: The evolution of the action density in a 1 degrees-of-freedom example. The red 
crosses are a quasi equilibrium while the green x's are a 0-dependent fluctuating density. The 
evolution is over 10 5 turns. 



5 The Linearized Equations 

To linearize about an equilibrium, we set ^(v,t) = \I/ e ( J") + *ffi(v,t), and ^/*(v,t) = ^l(J) + 
^i(v,t). Plugging into ( 1241) and linearizing gives 

/ + c{*i, h[%}} + a^e, Him = o, (9R) 

\ d&l + C{^1, H*[^ e }} + (*{%, H*[Vi]} = 0. 1 ' 

Because of the convolution structure of H in (1271) . the Fourier modes of (1251) are uncoupled, 
but we do not pursue fl28l) in that generality here. Instead we investigate solutions, which exist 
when m e = = (*, and H = H*. Letting = *i ± tfj, we obtain 

drf* + {/* H[^ e ]} ± {* e , = 0, (29) 

where f ± are the so-called a (sum) and n (difference) modes respectively and we have scaled 
the ( into the time by defining r = (t. We remind the reader that the densities and ^ 
in (128]) are coupled, while / + and /~ in (J29J) are not. Moreover, the a and 7r modes can 
be intuitively interpreted as in-phase and 180° out-of-phase perturbations respectively to the 
starred and unstarred beam. Note that t in (|28|) is dimensionless (same as n in (|22|) ) whereas 
t in ([2"9l has the dimension of a length. Note also that the beam current, which is a factor of 
(, has been scaled out of the problem. Thus there can be no threshold for instability in our 
problem. 

Unfolding the Poisson brackets we obtain 

d T f± + n(J) ■ V e / ± T Vj^ e ( J) • V e ^[/ ± ] = 0, (30) 

where 

n(J) :=VjH[V e }(J) . (31) 
To analyze equation fl30|) . we expand 

/ ± (J,e,r) = ^/ fe ± (J,r)e* e (32) 

fcez 2 
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To obtain an equation for the Fourier coefficients, we define 

K(J,J',e) := G (D(f3 x J X} f3 x J' x} Q x ), D(f3 y Jy, f3 y J y , <d y )) (33) 
with Fourier coefficients given by 

K ^ J ') = TTY2 [ K U J 'i Q)e' lhe dQ. (34) 

Since K is real and even in Q x and Q y , K k is real, and since K k is real, K k = K- k . Also 
K(J, J', Q) = K(J' : J, 6), and 

G(D X , D y ) = Kk(J, J') e ik < e - & \ (35) 

An easy calculation using (135]) and (l2Tj) gives 

H\f*\{J, 6) = (2vr) 2 / 2 K k (J, J')ft(J') df. (36) 

Thus the Fourier modes determined by (|29|) are uncoupled and are given by 

- id T fi + A; • £1( J)/f 
Tk-VUJ) [ K k (J,J')f±(J')df = 0, (37) 

where / e (J) := (27r) 2 \I/ e ( J) is the equilibrium action density ( L 2 f e (J)dJ = 1). 

There are two standard approaches to analyzing (1371) : the Laplace transform approach and 
the eigenvalue approach. 

Taking the Laplace transform of (I37p we obtain 

(co - k ■ n( J)) /±( J) 

±k-Vf e (J) [ K k {J,J')f±{J')dJ' = -ifjt(J,0), (38) 

for Su; sufficiently large. Here we use — icu instead of the usual Laplace variable s, f k (J,u) 
denotes the Laplace transform of f k , and Q(J) := H(^f e )(J)'. In the eigenvalue approach, we 
look for solutions of the form 

f k Ur) = f k ± (J)e-^\ (39) 

which gives 

(uj - k ■ n( j)) /±( j) 

±k-Vf e {J)l K k (J,J')f±(J')dJ' = 0. (40) 

The Laplace transformed equation fl38l) has a non-homogeneous term, however, the left hand 
sides of equations fl38l) and (j4TJl) are identical. These equations are of the form 



a(x)(j)(x) - X I K(x,y)(j)(y)dy = f(x) (41) 
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and are called Fredholm integral equations of the third kind (see p. 2 [16]). If a is bounded away 
from zero it can be transformed into a Fredholm integral equation of the second kind. Thus 
the primary interest in this third kind equation is in the case where a has zeros and this is our 
case. The case where a has zeros is complicated by the fact that there are generalized solutions 
which are difficult to represent numerically. 

Equations fl38l) and (HOl) have analogues in both plasma physics and other beam dynamics 
contexts. For example, Crawford and Hislop [T7] discuss the standard plasma problem in the 
periodic case, the case of this paper, summarizing both the Landau and the van Kampen-Case 
solutions ([18]. [19], [20] ). Jackson [21] gives a nice presentation of the Landau approach in the 
non-periodic case. The third kind integral equations are given in equations (17) and (23) of [17] 
and in equation (3.5) of [21] . As is well known, the plasma problem leads to Landau damping 
and growth for certain equilibria depending on the size of the average density. The stability 
analysis is facilitated by the dispersion function which uncouples the calculation of the poles of 
the solution from the calculation of the density itself. 

Two standard beam dynamics problems concern the longitudinal dynamics with wake fields 
for a coasting beam and for a bunched beam. The coasting beam case is completely analogous 
to the periodic plasma problem including a dispersion function and possible Landau damping 
for small beam current and an instability threshold at some critical current after which there is 
Landau growth. A recent discussion of the coasting beam problem in the context of coherent 
synchrotron radiation is given in |22j. The third kind integral equation is given in equation 
(27) of that paper. The emphasis in [22] is on the threshold for instability which occurs when 
a zero of the dispersion function reaches the real axis as the current increases from small 
values. Landau damping is not discussed as it is not important for the stability discussion and 
furthermore would require an analytic continuation into the lower half u plane which would 
require assumptions on the equilibrium (see p. 7 in [22]). 

The bunched beam case is more complicated as the Fourier modes do not decouple. Fur- 
thermore, it appears at first sight that the calculation of the instability threshold must be done 
in combination with the calculation of the density. However, Warnock has introduced a reg- 
ularization transformation which eliminates the continuous spectrum. The resulting equation 
is then discretized leading to a determinant condition, independent of the density, which is 
analogous to the dispersion relation. A convergence theorem would then make this rigorous. 
This is discussed in [23], where equation (11) is very similar to our equation (|38|) . However the 
kernel of the integral equation is much different and, in fact, one expects a stability threshold. 
It is in this paper that Warnock introduces his regularization transformation which eliminates 
the continuous spectrum and thus eliminates the numerical problem of representing generalized 
functions numerically. More recent progress on the regularization is given in [21]. 

Our equations fl38|) and fHUj) are simpler than the longitudinal bunched beam equations in 
that the Fourier modes are uncoupled. Also, our case is rather special in that it does not depend 
on the beam current as mentioned above. In fact the tt and a eigen-modes are neutrally stable 
if k ■ V/ e (J) 7^ 0. In this case, the transformation 



k-Vf e (J)\^ g±(J) 



(42) 



leads to 



fc.n(j))sf(j) 




(43) 
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where K k (J, J') := \k ■ V f e {J)\ l/2 K k {J, J')\k ■ V/ e (J')| 1/2 and tj* := ±sgn(£; • Vf e ). Since the 
kernel K k is real and symmetric, an eigenvalue u must be real. It is in fact a remarkable feature 
of the linearized AVE (I43p for the case of v = u*, [y x , u y ) far-from- low-order-resonance and for 
equilibria with k ■ V/ e (J) 7^ 0, that despite the presence of an amplitude dependent tune shift 
and a collective force, the modes show neither damping nor growth but instead are stable. In 
[261 [13] we have given numerical evidence that the modes are indeed remarkably stable even in 
the fully nonlinear regime of tracking with equation (J2J). 

We have tried to show that eigensolutions (to, f^) of PP|) for equilibria, which do not satisfy 
the above condition, k ■ V/ e (J) 7^ (e.g. densities with two humps), must have w 6 R, but 
have been unsuccessful. This leaves open the possibility of complex eigenvalues. Since these 
must come in complex conjugate pairs, there is the possibility of linearly (and thus nonlinearly) 
unstable solutions. 

We have assumed that the two beams have the same nonresonant tunes and this is probably 
the reason that the eigensolutions are neutrally stable. Previous work indicates that when the 
tunes are different (a so-called tune split) or near-to-low-order resonance there can be Landau 
damping or growth. In [25J, the authors develop a perturbation procedure in the near resonance 
case and argue that there are regions of stability and instability (see Figure 2 of [22] )■ Landau 
damping is also discussed in [1]. In [25], we have seen evidence for Landau damping in the 2D 
phase space case (See Figures 5-7 of[26j ). A future problem for us is to determine the explicit 
form of the averaged equations in this case. 

Equation (1431) and the associated Laplace transformed equation can be rewritten as (T — 
ujI)^ = h. Since T is symmetric, we are looking for conditions which ensure that it is a 
bounded selfadjoint operator on an appropriate Hilbert space. Such operators have a well 
developed spectral theory. For example, the spectrum is a compact subset of the real line 
contained in the interval [m, M] where both m and M are finite spectral values and all spectral 
values are either in the point spectrum (eigenvalues) or the continuous spectrum, thus the 
residual spectrum is empty [27]. Numerical results, in the section to follow, indicate that for 
k T = (1, 0) or (0, 1) the spectrum is the interval [0, k ■ fi(0)] for the a mode with an eigenvalue 
and [0, k ■ fi(0)] U {oj n } for the n mode, with u n > k ■ $7(0) an eigenvalue. A possible explanation 
for the stability of the modes using the notion of the "Landau resonance" is that the modes can 
not resonantly exchange energy with a macroscopic fraction of the particles in the beam. The 
a mode tune lies at the edge of the incoherent (continuous) spectrum towards infinite orbital 
amplitudes. Any sensible phase space density falls off rapidly at large amplitudes (or even has 
compact support) so that the fraction of particles with tunes in resonance with the a mode 
tune vanishes. The tt mode tune, in the studied parameter regime, lies considerably outside the 
incoherent spectrum and is thus even less able to dissipate its energy among the single particle 
trajectories or to draw energy from them. 

6 Numerical Results for tt and a Modes 

The analysis of the spectrum for equation (j43p gives important information about the it and 
a modes. In this section, we discuss properties of solutions of (1431) and give our results on the 
numerical solution of this eigenvalue problem. 

If to is outside of the range of k ■ fl(J), ff4"3"j) can be reduced to an integral equation of the 
second kind by a simple algebraic transformation. Conversely, if u is in the range of k ■ fl(J), 
then (j4"3"|) it must be treated as a third kind equation. Such equations have not been studied as 
extensively as Fredholm integral equations of the first and second kind. A review of work up 
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to 1973 and new results are given in [28J and more recent results are contained in [29]. Most 
recently, we have become aware of [3U1 ED E2] • However, to our knowledge, the case when J is 
2-dimensional has not been discussed nor have convergent numerical schemes been developed. 
As mentioned in the previous section the plasma problem and the coasting beam problems are 
of this type and have been studied extensively, however these are particularly simple. 

We now discuss the commonly used, straightforward discretization for integral equations of 
the third kind as applied to our special case and give our numerical results. At the end of this 
section we will discuss progress on work toward a convergent scheme. In the straightforward 
approach, J is put on a mesh and the integral is approximated by a simple quadrature method. 
This leads to a finite dimensional matrix eigenproblem and seems to lead to reasonable results. 
This approach has been used in ID in the beam-beam interaction in [31 H], in the longitudinal 
bunched beam case in [33] and by us in the ID beam-beam interaction, [13]. In [13] we obtained 
excellent agreement between the FFT spectra of the dipole modes in full blown simulations and 
the eigenvalues of a one degree-of-freedom version of (1431 . 

We consider the special case of axially symmetric Gaussian beams, where 

/e(J) = ie"^ (44) 

with ire being the rms emittance, and horizontal dipole modes where k = (1, 0) T . With these 
choices equation (l4"3]) takes the explicit form 

(uj - n x (J)) gf fi (J) 
Te- 3 [ e-^K^J'y-^gf^dJ' = 0. (45) 

We transform the actions I x = J x / (1 + J x ) and I y = J y / (1 + J y ), thereby mapping M + — > [0, 1), 
and use a 60x60 mesh. Even though the linearized averaged Vlasov equation ( 1451) reduces the 
number of independent variables from four as in ( 1241) to two, the evaluation of the functions 
Q( J) and Kk{J, J') is in fact computationally expensive. The computation of Q = VH[^ e ] 
involves a 6-fold integral at each point of the 2D mesh in J and involves a 4-fold integral at 
each point of a 4D mesh in (J, J'). Although we found a way to slightly simplify the calculation 
for general \l/ e , and reduce the 6- fold integral to a 5-fold integral, going to larger meshes is quite 
expensive. However, for the important particular choice of ( ]44"1) we found a very simple formula 
involving modified Bessel functions: 



dq 



2e + q 



2e + q J \2e + qJ \2e + q 

This formula has been known in the context of the weak-strong approximation of the beam- 
beam tune shift, [M]. It is straightforward to prove that limj^o^x(^) = hnij^o ^y(J) — 1> 
and that the ranges of Q x ( J), i\(J) are both the interval (0, 1). The latter is also the range of 
the continuous spectrum of (T43|) . 

^x{JxiJy) is shown in Figure [2j and the spectrum of the finite dimensional approximation of 
( |45l) is shown in Figure [3] The plot suggests that (]4"5l) has a continuous spectrum, common to 
both the a and 7r modes, which coincides with the range of Q x ( J). In addition, the a mode has 
a discrete eigenvalue u = 0, and the 7r mode has a discrete eigenvalue at u; ~ 1.21. Figures H] 
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Figure 2: Q X (J X , J y ) for the equilibrium density (jUj). 
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Figure 3: The dipole spectrum of fj45l) . The black curve is the continuum, common to both 
modes, the red circle is the discrete eigenvalue for the 7r mode and the blue diamond is the 
discrete eigenvalue for the a mode. 
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Figure 4: The a mode eigenfunction of ( |45l) . 
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Figure 5: The 7r mode eigenfunction of fj45|) . 
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Figure 6: An eigenfunction corresponding to an eigenvalue from the continuous spectrum of 
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and [5] show that these eigenvalues corresponds to regular eigenfunctions. FFT spectra obtained 
by tracking, with our newly written parallel Perron- Frobenius code which tracks the phase 
space densities directly in 4D phase space, and with our 4D weighted macro particle tracking 
code [IB], shows pronounced peaks at tunes that correspond to these discrete eigenvalues. This 
indicates an excellent agreement between these three completely different approaches. We 
consider this a major feat. 

Even though our results are quite satisfactory, we are interested in a convergent procedure. 
A well established theory (See e.g., [35]) guarantees that the above described finite dimensional 
approximation of the integral equation converges as long as the operator is compact. However, 
as we mentioned before, Figure [3] indicates the presence of a continuous spectrum, which a com- 
pact operator can not have. In addition, the numerically computed "eigenfunctions" associated 
with the continuous spectrum show a singular behavior, as illustrated in Figure El This is to be 
expected [281 EH] and thus these third kind integral equations are not well suited for numerical 
analysis in their standard form. Therefore further numerical and analytical analysis requires a 
special treatment of equation (143]) . Such difficulties are common for equation (T4T|) where a(x) 
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vanishes at least once inside of its domain. The recent work for the longitudinal bunched beam 
Vlasov equation in [231 Elj mentioned above, suggests that the continuous spectrum can be 
eliminated by the Warnock regularization transformation g(J, u) = (u — k ■ Q(J))g(J, to). This 
gives the nonlinear eigenvalue problem in u: 



where the solutions (g(J,u),u) for 77 = ±1 give the a and 7r eigenmodes respectively. 

As mentioned before, equation (T4T]) appears in other physics applications, and analytical 
results have been published in [28], [291 EH EH [32] . We have begun a study of this problem in our 
particular case. Specifically, the existence and uniqueness of solutions of equation (T4"3"l) in the 
one and two-dimensional cases has been addressed in [SB]- In [36], we consider functions which 
are continuous except for pole-like singularities at J such that (a? — k ■ £l(J)) = 0, and interpret 
the integrals over the singularities in the principal value sense. Under certain assumptions, we 
proved a version of the Fredholm alternative theorem: the equation has a unique solution for 
any right-hand side iff the homogeneous version of this equation has only the trivial solution. 
In addition, this framework may provide a numerically convergent scheme for solving (T4"3"j) . 
Alternatively, we are looking for conditions so that discretization of the Warnock transformed 
problem will lead to a convergent method. The theory of [36] is a first step. 

7 Summary and Future Work 

An important aspect of collective beam-beam theory has been the study of the so-called 7r and a 
modes. The pioneering works of [H [2], [3[ H] represent a major advance; however approximations 
are made, the validity of which we would like to understand better. For example, in [3j H], 
the starting point is a Vlasov equation with a delta function kick and the beam-beam kick is 
distributed around the ring by smoothing the delta function. The Vlasov equation is linearized 
about a function which is only an approximate equilibrium and action-angle variables are 
introduced. The Fourier modes in angle are not uncoupled at this stage but a horizontal dipole 
mode proportional to exp(iQ x ) is assumed. This leads to an inconsistent equation. To obtain 
a consistent equation for this mode an average over the angle variables is taken which, after 
a Fourier transform in time, leads to an integral equation, the analog of equation (T4"3"l) . In 
contrast, we start with the kick-lattice model to properly handle the delta function kick. Then 
we make only one approximation, the averaging approximation, and in addition, as stated 
earlier, we believe we can give an upper bound on the error of approximation. Our AVE has 
exact equilibria and thus our exact problem has quasi-equilibria in good agreement with our 
simulations. The linearization about these equilibria leads to an equation, which in contrast 
to the above, has uncoupled Fourier modes. The Fourier modes satisfy an integral equation 
that is easily transformed to a formally selfadjoint problem. We are looking for conditions such 
that the associated operator is bounded and selfadjoint, a case which has a well developed 
theory. The standard computation of the 7r and a dipole mode frequencies discussed in Section 
6 is in good agreement with density tracking based on equation (J3j), using both the PF and 
the weighted macro-particle tracking methods. However the standard numerical approach to 
numerically solving (1431) does not converge as the mesh size decreases beyond some limit and 
we are searching for convergent algorithms such as that suggested in [23l [36] . 

In summary, we have introduced a new model, the averaged Vlasov equation (1241) . for the 




(47) 
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collective beam-beam interaction in two degrees of freedom, which we believe has significant 
potential for deepening our understanding of this important collective effect. 

Equation (j24j) was derived in the spirit of the rigorous analysis in [H] . We believe similar error 
bounds can be derived, thus we believe (jUj) gives a good approximation to the basic dynamics 
of (j2J). In fact, we have checked the one degree-of-freedom analog of ff24l) with two aspects of 
a full-blown density tracking approach, the existence of quasi-equilibria and the calculation of 
the it and a mode eigenvalues, with excellent results. More importantly, we have checked the 
two degree-of-freedom AVE with two full-blown simulation codes and have also found excellent 
agreement in the calculation of 7r and a mode eigenvalues. Thus we have confidence in the 
model. 

We have demonstrated its usefulness as a tool for calculating 7r and a mode eigenvalues 
and for clarifying the existence of quasi-equilibria. In the case of leptons, progress has been 
made on the question of the existence of an equilibrium for the exact model [T5]. However, it 
seems likely that exact equilibria do not exist in the hadron case as the underlying dynamics 
is likely to be chaotic. In addition, the AVE may lead to a faster algorithm for calculating the 
density evolution. This is because the beam-beam parameters, ( and (*, are small and thus the 
time step in numerical integration of (EMI) can be 0(l/max(C, C*)) ; which is significantly larger 
than one turn. We propose to investigate this potential speed up by developing a numerical 
procedure to integrate (I!MI) . As another example, the AVE will be useful in taking the next 
step beyond the linear theory to investigate coupling between the n and a Fourier modes. In 
the Laplace-picture, we may be able to use the fixed point iteration scheme discussed for the 
plasma problem in [37] or in the eigen-picture presented here we may be able to use the van 
Kampen-Case approach [T9l [20] . In the latter case, the work of [17J may be useful. Finally, we 
can investigate several other effects such as those discussed by Alexahin [I]. Some topics we are 
considering for future work are: (i) a study of the near-to-low-order resonance case as we do 
in [6] (ii) adding another degree of freedom to study the effect of synchrotron motion into the 
dynamics of equation (2) and (iii) a study of the effect of a tune split by letting u* — v x + (a 
and u* — v y + (b and then applying our averaging formalism. Here (a, b) allow us to vary the 
tune split in units of the kick parameter. As in [5] we expect bifurcations as a and b vary. Items 
(i) and (iii) likely includes the possibility of both Landau damping and growth. 

Our main point is that we now have a model in which many important collective beam-beam 
interaction effects can be studied in a more systematic way then was previously available. 
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A Derivation of the Beam-Beam Parameter 

To calculate the kick we consider three inertial reference frames: the rest frame of the syn- 
chronous particle of the unstarred bunch F, the lab frame F', and the rest frame of the syn- 
chronous particle of the starred bunch F*. The coordinate axes of the three frames are parallel 
and oriented so that viewed from the lab frame F', F is moving in the positive z direction with 
velocity (3c and F* is moving in the negative z direction with velocity —f3*c. In what follows 
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we will consider the relative velocities of the particles with respect to the synchronous particle 
of their bunch as non-relativistic. 

We consider the momentum change of an unstarred particle moving through the starred 
bunch. In the kick approximation, we assume that the transverse spatial coordinates are not 
changed during the interaction, i.e. r(t) = (x,y,z + ut) T . Thus the change in transverse 
momentum of an unstarred particle passing through the starred beam measured in F* is 

APl(x,y) = q [ E* ± (x,y,z + ut)dt=l [ E* ± (x,y,z)dz, (48) 
Jr u Jr 

where u is the speed of the unstarred particle in the starred frame, 

u = c — — , 49 

and is the transverse electric field of the starred bunch in F*. Note that j3 is basically the 
speed of an unstarred particle in the lab frame and that B* is approximately zero in F*. 

Since the 3-momentum P* is part of a 4-vector and the boosts involved are all in the 
longitudinal direction 

AP ± = AP f ± = API. (50) 

From = — V_l0* we have 

APl(x,y) = -^*- [ dzV_ 

xpl(f ± ,z)-== (51) 



N*qq* 



/ P* 3 (r±,z) / dzVj 
Jr 3 Jr 



x - £M (52) 
V \\r± - r-_L || 2 + (z - zf 

One easily shows by direct integration that 

dzV ± ) = -V ± In ||r ± -f ± f, (53) 

yf\\r ± - r ± \\ 2 + {z- z) 1 

which is independent of z. Thus 

APl(x, 2 /) = ^^Vx / P ;(rs_)ln\\r ± -r ± \\, (54) 
uAneo J K 2 

S v ' 

:=<J>[p*](r ± ) 

where p\ = L p% dz. The interchange of limits going from fl5T|) to (1521 is justified for p% decaying 
sufficiently fast at oo. (The singularities in the integral are integrable.) 
Since AP'j_ = AP± the kick in the lab frame is 

A P y J p Q 
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Since the speed of the unstarred synchronous particle in the lab frame is (3c, we have po = m(3^c. 
Thus the kick is 



Ap. 



11 



-CV ± m\(r±) (56) 



where 

= 2N*qq* J_ l + (3(3* 
^ 4nme c 2 /3 7 (3 + (3* ' 1 j 

which is the justification for equation ([2]). 
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